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ABSTRACT 

Radiative transfer calculations of massive star formation are presented. These are based on the 
Turbulent Core Model of McKee & Tan and self-consistently included a hydrostatic core, an inside-out 
expansion wave, a zone of free-falling rotating collapse, wide-angle dust-free outflow cavities, an active 
accretion disk, and a massive protostar. For the first time for such models, an optically thick inner 
gas disk extends inside the dust destruction front. This is important to conserve the accretion energy 
naturally and for its shielding effect on the outer region of the disk and envelope. The simulation 
of radiation transfer is performed with the Monte Carlo code of Whitney, yielding spectral energy 
distributions (SEDs) for the model series, from the simplest spherical model to the fiducial one, with 
the above components each added step-by-step. Images are also presented in different wavebands of 
various telescope cameras, including Spitzer IRAC and MIPS, SOFIA FORCAST and Herschel PACS 
and SPIRE. The existence of the optically thick inner disk produces higher optical wavelength fluxes 
but reduces near- and mid-IR emission. The presence of outflow cavities, the inclination angle to 
the line of sight, and the thickness of the disk all affect the SEDs and images significantly. For the 
high mass surface density cores considered here, the mid-IR emission can be dominated by the outflow 
cavity walls, as has been suggested by De Buizer. The effect of varying the pressure of the environment 
bounding the surface of the massive core is also studied. With lower surface pressures, the core is 
larger, has lower extinction and accretion rates, and the observed mid-IR flux from the disk can then 
be relatively high even though the accretion luminosity is lower. In this case the silicate absorption 
feature becomes prominent, in contrast to higher density cores forming under higher pressures. 
Subject headings: ISM: clouds, dust, extinction — stars: formation 



1. INTRODUCTION 

"How do massive stars form? " is still a debated ques- 
tion (e.g., Beuther et al.|[2QQ7l ). One basic problem is 
massive protostars become so luminous that radiation 
pressure may stop the accretion and growth of the star. 
One possible mechanism to form a massive star can be 
considered as a scaled-up version of low-mass star for- 
mation. If the gas is turbulent and threaded by suffi- 
ciently strong magnetic fields then fragmentation may be 
suppressed for cores much more massive than the Jeans 
mass. The conditions leading to this suppression should 
be relatively rare since massive stars are rare and make 
up only a small mass fraction of the final star cluster. 
These massive cores are expected to form from highly 
pressurized clumps of gas, in which case they start with 
high densities, short free-fall times and therefore high 
accretion rates. This is the basic scen ario of the Turbu- 
lent Core Model (jMcKee fc Tanll2QQ3L hereafter MT03). 
Other radically different possibilities are tha t massive 
stars form through stellar mergers (B onnel etld]ll998f ) 
or accrete most of their mass from init ially unbound ma- 
terial (competitive accretion n iodel, iBonnel et al.ll2QQll 
iBonnel et £0112004 lBatell2QQ9al lbh. 

Answering this question is difficult observationally be- 
cause massive star formation occurs in distant and highly 
obscured regions. However, with the improving sensitiv- 



ity and spatial resolution provided by the instruments 
such as the Herschel Space Telescope^ Stratospheric Ob- 
servatory for Infrared Astronomy (SOFIA) ^ Gran Tele- 
scopio Canarias (GTC) - CanariCam^ the James Webb 
Space Telescope and Thirty-Meter class telescopes, we 
expect to see a faster advance in the research of mas- 
sive star formation and hope this question can be finally 
solved. 

In order to interpret observations such as spectral en- 
ergy distributions (SEDs) or images, and then have a 
better understanding of the properties of a massive pro- 
tostar and its evolution, a number of models have al- 
ready been develop ed to fit or coni p are w ith observa- 
tions. For example, iRobitaille et al.l (|2006f ) have devel- 
oped a very impressive model grid containing 200,000 
SEDs covering a large parameter space, which is pub- 
licly available and now widely used (referred to below 
as the "Robitaille model"). However, this grid mainly 
covers low-mass young stellar objects (YSOs) and when 
massive protostars are considered, they do not have the 
properties expected for fiducial parameters of the tur- 
bulent core model. Also these models do not consider 
the presence of an opticall y thick gaseous disk inside 
the dust destruction front. iMolinari et al.l (|2008l ) used 
the same methods employed by the Robitaille model but 
now for protostellar parameters based on the turbulent 
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core model, to developed a SED model grid for massive 
YSOs. They found that the SED can be a diagnostic 
tool to determine the evolutionary stage of a massive 
YSO. However, in this work the representation of the 
turbulent core model is quite approximate and limita- 
tions of the Robitaille niodel framework are still present. 
IChakrabarti fc McKed (|2QQ5f ) developed an analytic so- 
lution for far-IR SED of a protostar embedded in a spher- 
ically symmetric molecular cloud, but this method could 
not allow for presence of a ccretion d isks and protostellar 
outflow cavities. Indebe touw et all (12006) investigated 
the effects of clumpy structures in the molecular enve- 
lope around massive protostars. However, again, their 
models were not tuned to the parameters of the tur- 
bulent core model: for example, they considered much 
more massive structures (e.g. ~ 5 x 10^ Mq contained 
in a sphere of radius 2.5 pc), more representative of star 
forming clumps that form entire star clusters. They did 
not consider protostellar disks. 

Our aim is thus to develop a new model of mas- 
sive protostars , base d on the Turbulent Core Model of 
iMcKee TanI (|2002[ ) and MT03, including all the im- 
portant components self-consistently, and then perform 
simulations of the radiation transfer to see whether dif- 
ferent components or evolutionary stages are represented 
in the SEDs and images. Starting with a fiducial hydro- 
static core of 60 Mq bounded by the pressure of a self- 
gravitating clump of mean mass surface density Hci = 
1 g cm~^, we then consider its appearance once an 8 Mq 
protostar has formed at its center. We develop a se- 
ries of protostellar models of increasing realism: starting 
with a simple hydrostatic cor e, we the n apply the inside- 



out expansion wave solution (|Shul 19771), but generalized 
to si ngular polytropic spheres (|McLaughlin fc Pudritz 



|1997) and the free-fall rotating collapse solution bv lUlrich 
fl976). A circumstellar disk is expected to form around 
the protostar and this is important to transfer angu- 
lar mome ntum and to solve the r adiation pressure prob - 
lem (e.g. iJiiina fc Adanislll996L fKrumholz et"ani2007[ ). 
Presently, there is evidence for rotating toroids around 
massive protostars (e.g. Beltran et al. 2005) but little di- 
rect evidence for Keplerian protostellar disks, which will 
probably require the angular resolut ion of ALMA. We in- 



clude the disk with an a-disk m odel fShakura fc Sunvaev 
Unlike iRobitaille et al.l (1200 6*) and' Molinari et al 



(|2008f ). we include an optically thick inner disk with 
gas opacities inside the dust destruction radius, which 
is important to conserve the accretion energy naturally. 
Accretion is expected to drive strong bipolar outflows 
and sweep up the material in the core and form cavi- 
ties. These outflow cavities have been observed around 
massive protostars and may determine the mid-IR mor- 
phology (e.g., De Buizer 2006) . 

The assumptions of each component of our model and 
the model series are introduced in detail in the next sec- 
tion. In Section. [3l we discuss our simulations, including 
the Monte Carlo radiation transfer code, and the dust 
and gas opacities we use. In Section. IH we present the 
SED and image results of our models. In Section. O we 
summarize our main results, including a comparison of 
our model with other works. In future papers we will ex- 
amine additional refinements, especially the development 
and material content of the outflow cavities, and present 
results of the evolutionary sequence of massive star for- 



mation based on the fiducial model we have started to 
develop here. 

2. MASSIVE PROTOSTAR MODEL 

2.1. Envelope 

2.1.1. Hydrostatic Outer Envelope 

Following MT03 and (|2008f ). we define a "star- 
forming core" as a region of a molecular cloud that will 
form a single star or close binary, and assume it is self- 
similar, self-gravitating in near virial equilibrium and 
spherical. The density and pressure each have power-law 
dependencies on radius, p ex r~^p and P oc r~^p . This 
smooth power-law density distribution is only an approx- 
imation, especially given that massive cores in virial equi- 
librium but with gas temperatures ~ 10 K must be sup- 
ported by non-thermal forms of pressure support, such 
as turbulence and/or magnetic fields. Clumpy substruc- 
tures are likely to form inside a turbulent core and this 
will affect radiative transfer through the core - basically 
making it somewhat easier for shorter wavelength pho- 
tons to propagate through the core. Ilndebetouw et afl 
(i2006i) performed radiative transfer models of clumpy 
cores. However, given the uncertain nature of the clump- 
ing, we defer such considerations to a future paper, and 
first calculate the properties of smoothly distributed gas 
and dust. 

From the above power law distributions, it follows that 

the core is polytropic with P oc p^p. T he ca se with 

7p = 1 is a singular isothermal sphere (|Shulll977l ) and in 
the other limit 7p = corresponds t o a logotropic sphere 
(McLaughl in fc Pud ritz"1996. 1997). For the models pre- 
sented here, we follow MT03 and adopt kp = 1.5, thus 
7p = I and kp = 1. The equation of hydrostatic equilib- 
rium gives 



M(r) 



kpC^r 



and 



p{r) 



(3- 



l)c2 



27rGr2 



(1) 



(2) 



where c = (P/p)^/^ is the effective sound speed. 

We assume the total mass in the core is 60 Mq. If 
the efficienc y e^f is 0.5, which is es timated from low- 
mass cases (Matzner fc^lcKed"2000), a star with mass 
m*/ = e*/Mcore = 30 Mq can finally form out of the 
core. For a core with such mass, the core radius is (eq. 
(20) in MT03) 



(3) 



where is the mean surface density in the molecular 
clump in which the core is embedded, and 1 g/cm~^ is 
used as a fiducial value. We also consider models with 
higher and lower values of Sci- 
2.1.2. Expansion Wave 

IShul (|l977l ) developed the inside-out expansion-wave 
solution for the probl em of the gravitational collapse of 
an isothermal sphere. iMcLaughlin fc Pudritzl ([l997) ap- 
plied this solution to the collapse of a logotropic sphere 
and also gave the general formulae for a polytropic 
sphere. Following their work, we calculate the density 
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Fig. 1. — Density distribution with radius r (left panel) and polar angle (right panel). riHe = O.lnn is assumed here. In the left 
panel, different curves correspond to different polar angles, which is 0° when along the rotation axis. The thick black curve is without any 
rotation effect for comparison. The right panel shows the density dependence on ^ at r = 2rc/ and r = Ard. 



profile in the core. A core with initial density profile of 
eq. ([2]) is not stable and a perturbation at the center 
can trigger collapse of the innermost region. The posi- 
tion where the material begins to fall progresses outward, 
which is called the expansion wave. While the material 
collapses inside the expansion wave front, the region out- 
side it is still hydrostatic. The similarity is always kept 
during the whole collapse in this solution. 

iMcLaughlin fc Pudritzl (|l997l ) did not include the ef- 
fect of magnetic fields which can help to support more 
mass in the core and increase the accretion rate after 
collapse starts. We estimate the effect of magnetic fields 
from the work of Li & Shu (1997), who found that the 
equilibrium surface density is increased by a factor of 
(1+i^o) when magnetic fields are considered, where Hq is 
a parameter. So the real mass and density profile should 
increased by the same factor, 

M(r,t) = Mnon(r,t)(l + i7o) = m{x){l^ Ho)a^tt/G, (4) 

and 

-r^ N a(x)(l -\- Hq) . , 

p{r.t) = P„o„(r,t)(l + ffo) = ^ ^^^^^2 ' (5) 

where m{x) and are the similarity variables for mass 
and density, and at = [i^7p(47rGt^)^~^p]^/^ has dimen- 
sion of velocity. In our case, i^o is set to 1 following the 
assumption by MT03. 
The collapsed mass at the center now is 

M(0) =m(0)(l + i^o)a?t/G, 



therefore, 
and 



(6) 
(7) 
(8) 

in our case, which is consistent with the analysis of 
MT03. So, given a certain collapsed mass M(0,t), with 
the star formation time (the time that the whole core 
takes to collapse, eq. (44) in MT03) 



M(0) (X a^t (X 
M(0) oc t^-^^^ (X t 



and the final collapsed mass M(0,t = t*/) = Mcore = 
60 M0 (The whole core collapses at the end, either into 
the star-disk system or into the outflow and escapes), we 
can calculate t and further the density p{r^t\ velocity 
u{r^t) and mass M{r^t) at that moment. At some time 
the expansion wave will reach the boundary of the core 
and lead to a backward wave, thus, for simplicity we 
choose such a collapsed mass that the expansion wave 
has not reached the core boundary yet. (We will discuss 
this in detail in Section [2^ ) 

2.1.3. Rotating Infall 

We consider a slowly rotating core. For simplicity, we 
only include the effect of rotation inside the sonic point, 
where the infall becomes supersonic. Once such infall 
starts, we expect it is difficult to transfer angu lar momen- 
tum f rom inside the sonic point to outside (Tan fc Mc Ked 
[2QQ1. So we assume that the angular momentum is con- 
served inside the sonic point until gas accretes onto the 
star or disk. 

We use the solution of lUlrichl (|l976f ) (referred to below 
as the "Ulrich solution") to describe the velocity field 
and density profile inside the sonic point. This solution 
assumes that a particle with an initial distance Too from 
the center, an initial polar angle ^0 and angular velocity 
Vtoo about the axis of rotation, due to a point mass M 
at the center, moves in a parabolic path described by 



Vd cos c/o sm 6/0 
cos ^0 — cos 



where 



(10) 



(11) 



(9) 



GM CM 

All particles starting from a spherical shell of radius Too 
will hit the equatorial plane at radius r < r^i, so naturally 
Vd can be thought of as the radius of the accretion disk. 

Td is also the centrifugal radius which marks the radial 
extent of centrifugal b alance, inside which t he flow will 
become disk- like (e.g. iJiiina fc Adamslfl996[ ). It can be 
estimated as below: 

ra='^Pr='-^p(-^)^, (12) 
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where ag = 5|^grav|r/(3GM2) = (5/3)(3 - kp)/{b - 
2kp) 1.25, ai = 2E,ot/{Mr'^Q?) = (2/3)(3 - /cp)/(5 - 
kp) — > 0.286, and /3 = £^rot / 1 ^grav | is set to have a fidu- 
cial value of 0.02 for the slowly rotating core, based on 
observations of low-mass cores (jGoodman et al.l 11993). 
For any spherical shell at radius r inside the sonic point, 
it infalls following eq. (fTQ|) , with rd to be the outer radius 
of the disk determined by the position of the sonic point. 



6ag 



/ Msp \ 3-kp 



(13) 



where Mgp is the mass inside the sonic point. In our 
fiducial case, when the collapsed mass is 10.67 Mq (8 
Mq star and 2.67 Mq disk, discussed in Section [23]), 
the disk radius is 449 AU. 

The rotation changes the density distribution from 
spherical symmetry to axissymmetry as below 



cos^ \ -1/2 / cos^ 



, / cos(y\-v^/ 

p{r,0) (X 1 + 

V cos Oq/ \ 



^ 2 cos 00 



^ COS^ ^0 

r 



(14) 

This result is valid for an infall rate that is constant with 
radius, which is not exactly right for the expansion wave 
solution of a polytropic sphere. Thus we only scale the 
density distribution with the angular dependence of eq. 
(p!4j) and keep the total mass in each shell unchanged, so 
that the infall rate is the same as the non-rotating case. 
Fig. [T] shows the radial density profiles at different polar 
angles (upper panel) and the density profiles with po- 
lar angle at different radii (lower panel), for a core with 
10.67 Mq collapsed at the center. In the upper panel, 
the thick black line shows the density profile in a non- 
rotating core. The thin curves correspond to different 
polar angles. There is a discontinuity at r ~ 2500 AU, 
which marks the position of the sonic point, only inside 
of which we consider rotation. For a real core, this dis- 
continuity would not exist. The peaks and dips mark 
the disk radius Vd. The lower panel of Fig. [T] shows the 
dependence of density on 6 at 2rd and Ard. 

2.2. Accretion Disk 

Inside the centrifugal radius rd, the infalling fiow is 
circularized and forms an accretion disk around the pro- 
tost ar. The disk can provide an efficient way to transfer 
angular momentum outwards and enable high accretion 
rates to form a massive star. The high accretion rate in- 
dicates that the disk is relatively massive (comparable to 
the star). But gravitational instability can become very 
efficient when the mass of the disk is high enough and 
lead it to fragment. We assume the disk mass is always 
a constant fraction of the stellar mass (|Tan fc McKed 
[200I : 



m^d = m* + = (1 + /d)m*. 



(15) 



where fd is assumed to have a fiducial value of 1/3. Re- 
cently, Kratter et al. (2010) performed a numerical pa- 
rameter study on accretion disks of massive stellar sys- 
tems and found that a disk can be stable and does not 
fragment with an even higher disk-to-star mass ratio 

{fd ~ 1). 

We follow th e description for the disk structure by 
IWhitnev et al.l (l2003b), which is a standard a-disk 
(jShakura fc Sunvaevlll97l iLvnden-Beh fc Pringlelll974 



iPrindd ll98lL iBiorkmanI I1997L iHartmann et al.l I1998D 
with density distribution 



' = Po (1 



where r is the radial coordinate in the disk midplane, z 
is the distance from the midplane, and is the stellar 
radius. Basically the density has a power-law profile with 
radius (smoothed at the innermost region) and a Gaus- 
sian profile with the height. H is the scale height of the 
disk, 

H = H.{^y. (17) 

Thus, the disk structure is specified by following param- 
eters: disk mass rrid. disk inner radius rd in 7 

disk outer 

radius rd (the centrifugal radius), disk scale height Hq at 
the stellar surface, and the power indices a and p. 

If the dust is the only source of opacity in the disk, 
then it equals to a disk truncated at dust destruction 
radius rgub, inside of which the temperature is too high 
for dust to exist. However, in reality, the gas opacity at 
the innermost region of the disk can be significant, so the 
disk should be optically thick down to the surface of the 
protostar. Thus, in our fiducial model, the inner radius 
of the disk is set to be the stellar radius. 

We consider both geometrically thin {Ho/r^ = 0.01) 
and moderately thick (Ho/r^ =0.1) disks to see the ef- 
fects of the dis k scale heig:ht. H o/r^ = 0.01 is the value 
used by (Whi tnev et al.ll2003al lbh. i^o/r* = 0.1 is closer 
to the aspect ratio expected for a disk composed of dust 
only and is a typical value for a moderately thick disk. 
For these cases, a and /3 are chosen to be 1.875 and 1.125, 
respectively, as in the standard a disk models considered 
by Whitney. However, in our final fiducial model, we 
calculate the disk scale height, a and P self-consistently 
from the gas and dust opacities, in which case, fitting 
with eq. (p!6|) and (pT|) gives Ho/r^ = 0.06 and a = 1.75, 
f3 = 1.08. 

In the limit of a slowly rotating protostar, the accretion 
luminosity from the system is 



Gm^m 



(18) 



half of which is emitted from the viscous disk (disk accre- 
tion luminosity Ldisk) and the other half is emitted when 
the material hits the surface of the protostar (hot-spot 
luminosity Lhotspot), 



-^disk — L 



1 



hotspot 



(19) 



For simplicity, we add the hot-spot luminosity to the 
star's homogeneously. The star then has a single, en- 
hanced temperature to describe its black body spectrum. 
If the disk is truncated at rd in, the disk luminosity will 
be 

W(r.,in) = ^^(3-2(r*/r,,in)'/'), (20) 

^'^(i,in ^ ^ 

which indicates that only when the disk extends to the 
stellar radius can the total energy be conserved. In fact, 
the disk accretion luminosity is very sensitive to the inner 
radius of the disk. The accretion rate m of a massive pro- 
tostar can be very high, reaching ~ 10~^ — 10~^ Mq/jt. 



Radiation Transfer of Models of Massive Star Formation. 1. 



5 



Here we adopt the protostar evolution model by MT03, 
which gives an accretion rate of 2.40 x 10~^ Mq/jt for 
an 8 Mq protostar in a 60 Mq core. This accretion rate 
gives a disk luminosity and a hotspot luminosity which 
are comparable to the stellar luminosity. 

The disk accretion rate is related to the a parameter 
of the disk by 

rfi = VlSn^a^iskVcPoHl/n (21) 

with Vc = (Gm*/r*)^/^. Parameters used in our fiducial 
model correspond to a case with adisk = 1.43, which is 
consistent with t he re sults of numerical simulations by 
iKrumholz et al.l (|2007[ ). in which they found an effective 
Q^disk ^ 1.0 - 1.6. 

Also, part of the energy may be used to drive the out- 
fiow (mechanical luminosity L^, MT03). But for now we 
only consider the radiative luminosity from the disk and 
the hot-spot. 

2.3. Outflow Cavity 

Powerful bipolar outfio ws are ubiquitous phen omena 
around protostars (e.g. Konigl fc Pudritzl l2000f ). The 
prevalent interpretation is that outfiows are powered by 
accretion activity, being driven by spinning magnetic 
fields that thread the disk. There are several theoreti- 
cal models to describe this process suc h as the X-wind 
from the innermost region of the disk fSh uet aLllTool 
and disk wind model (Konigl & Pudritz 2000f ). 

A common feature of these models is the produc- 
tion of a bipolar outfiow with momentum distribution 
oc (sin^^y)"^ for 6^ > O^Oi where 0^ is measured 
from the outfiow axis and O^o ~ 10 ~^ is a small angle 
([Matzner fc McKee.,1999) . On scales large compared to 
the source, 

dn ~ 4^1n(2/^^o)(l + - cos^ e^ ) ' ^ ^ 
We follow the discussion of iTan fc McKed (|2002f ). As- 
suming an opening angle of O^^qsc for the outfiow cavity, 
the outfiow material with 0^ > Ow,esc cannot escape from 
the core and will go back to the infalling fiow. We pa- 
rameterize the fraction of the out flow momentu m that 
escapes from the core with fw, esc- iNaiita fc Shul (p.994) 
showed that the velocity of the wind is approximately 
independent of the polar angle so that fw, esc can also 
describe the ratio between the mass loss rates. We as- 
sume that part of the material accreted to the disk will 
be transferred to the star at a rate m*, and another part 
will leave to the outflow at a rate m^^, while the rest is 
left in the disk so the disk grows in a rate rhd. We also 
assume that rhd = fd'^^i fw = rhw/rh^ and fw, esc are all 
constant. 

Generally, we assume the outflow starts at time ti 
when the stellar mass is m*(ti) and the disk mass is 
^d(^i) (^1 = corresponds to situation where the out- 
flow cavity forms as soon as the collapse begins.) When 
t < ti there is no outflow and we have 

^*d,o(^i) rn^dih) m^{ti)^md{ti) (l + /d)m*(ti), 

(23) 

where m^d,o is the collapsed mass of the polytropic core 
(a hypothetical star-disk mass if the feedback is absent). 
After ti until the time when the whole core has collapsed 
t/, we have 



where cos^^^^esc nieans only part of the infalling core ac- 
cretes to the disk because of the existence of the outflow 
cavity. Therefore, the instantaneous star formation effi- 
ciency is 

rn^d rn^djtf) - m^dih) _ (1 + fd)cosOw,esc 



- ^*d,o(^l) 1 + /d + fwfw ,esc 

(25) 

and the mean star formation efficiency is 



m^d(tf) rn^d(tf) 



i*d,o{tf) 



(26) 



After the whole core has collapsed, part of material in 
the disk will accrete onto the star and rest of them will 
leave the star-disk system to the outffow wind, which 



gives us 



(27) 

where m*/ is the mass of the star flnally born out of the 
core, which is assumed to be half of the initial total mass 
of the core Mcore, i-e, the flnal star formation efficiency 
is 



m^f fd^ fw 

Mcore " (l + /^)(l+/d) 



(28) 



+ + fw,e 



nTriw = COS( 



(24) 



For the models with outflow cavities that we investi- 
gate in this paper, we make the assumption that the out- 
flow cavity has only now formed when m*(ti) = 8 M©, 
so m^d{ti) = 10.67 Mq. The validity of this assumption 
will be examined in a future paper. Combining eq. (j22j) 
to ([28|) we solve for the opening angle O^^esc and fw, esc 
simultaneously. For e*/ = 0.5 and fd = 1/3, we flnd 
that when fw = 0.1 ~ 0.8, the opening angle varies from 
64° to 45° and /^,esc varies from 0.91 to 0.84. So the 
fraction of the outflow coming back to the infalling core 
is always small, and the opening angle is typically large 
except when fw is close to 1. We choose fw = 0.6 as 
our flducial value, making fw, esc = 0.86 and the opening 
angle of the outflow cavity to be 51°. In our models, 
the cavity wall follows the streamline of the rotating in- 
falling material (eq. [10]). For now we set the outflow 
cavity to be empty (i.e. the optically thin limiting case). 
One would expect it to be relatively free of dust if most 
of the outflow is launched from the region of the disk 
inside the dust destruction front and assuming there is 
little time for new dust grains to form in the rapidly ex- 
panding outflow. We will improve upon this assumption 
by studying the detailed density distribution in the out- 
flow cavity in our next paper. Since the cavity does not 
change the density distribution in rest of envelope, the 
total mass of the material in the envelope after outflow 
cavities are carved out is about half its original value. 

2.4. Protostar 

Assuming that the core mass is a constant, the mass of 
the central protostar indicates its age and evolutionary 
stage. At some moment, the outgoing expansion wave 
front will reach the boundary of the core and, depend- 
ing on the properties of the boundary, may lead to a 
backward wave and a breakdown of self-similarity. In 
our case, this happens when the collapsed mass reaches 
^*d,o ^ 12 Mq. For the present paper we have chosen 
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Fig. 2. — Structure of the massive star forming core in the fiducial model. The right panel is a zoom-in view of the central region of the 
left panel. The positions of core boundary, expansion wave front, sonic point (inside which we consider rotation), disk scale height, dust 
destruction radius and the outer radius of the disk are marked. The outer boundary of the disk is chosen to be a surface with constant 
density joining the outflow cavity wall at r^. 



to consider a series of models with m* = 8 M©, and thus 
a maximum central collapsed mass of 10.67 Mq for those 
cases with rotating infall to a disk. Thus the central ob- 
ject is on the verge of becoming a ^^massive protos tar" , 
following the definitions of Beuthe r et al.l (|2QQ7[ ) and lTanI 
|2008). The fuh evolutionary sequence, including a treat- 
ment for greater masses, will be considered in a future 
paper. 

MT03 studied the evolution of massive protostars, and 
for = 8 Mq and Ed = 1 g cm~^, the protostar is not 
expected to have yet contracted to the main sequence. 
They calculated the values of radius and luminosity of the 
protostar: r* = 12.05 Rq and = 2.81 x 10^ Lq. For 
simplicity, we assume a black body spectrum for the star, 
with surface temperature = 1.22 x 10^ K in this con- 
dition. For models with an active disk, we also add the 
hot-spot luminosity to the stellar spectrum, assuming it 
is emitted homogeneously from the stellar surface, which 
means the temperature now is T^^hotspot = 1.43 x 10^ K. 

For the cases with Hd = 0.316 g cm~^, we use the 
following parameters: = 11.3 i?©, = 1.25 x 10^ 
K, and T^^hotspot = 1.37 x 10^ K. For the cases with Sci 
= 3.16 g cm~^, these parameters are: r* = 5.93 i?©, 
T* = 1.74 X 10"^ K, and T*,hotspot = 2.62 x 10"^ K. 

2.5. Model Series 

In order to study the effects of all the features discussed 
above on SEDs and images, we construct a model series 
starting from the simplest to the one containing all these 
features. All the models assume an initial core mass of 
60 Mq from which an 8 Mq protostar has formed at the 
center. 

In Model 1, we assume a spherical symmetric density 
distribution in the core with a power-law dependence on 
the radius, p oc r^p . The total mass in the envelope is 
52 Mq since an 8Mq protostar has already formed at 
the center. 

In Model 2, We change the radial density distribution 



in the envelope to the expansion wave solution (similar 
to the thick line in the upper panel of Fig. [H but here 
the expansion wave front is at r = 0.0408 pc and the 
envelope mass is again fixed at 52 Mq). 

We begin to consider rotating infall inside the sonic 
point and thus a disk around the star in Model 3. The 
disk is geometrically thin {Hq/t^ = 0.01) and passive 
(no accretion luminosity) with the inner radius set to 
be the dust destru ction radius rsnb: wh ich is empirically 
determined to be (|Whitnev et al.|[2004f ) 

rsub =r*(Tsub/T*)"'-' (29) 
where Tgub is the dust sublimation temperature and we 
adopt Tsub = 1600 K. The expansion wave front now 
reaches r = 0.0494 pc for the collapsed mass now is m* + 
rrid = 10.67 Mq. The envelope mass is now 49.33 Mq. 

Outfiow cavities are added in from Model 4. We keep 
the density profile in the envelope unchanged, which cor- 
responds to a case that the outfiow has just swept up the 
material to form the bipolar cavities. The envelope mass 
is now ~ 29 Mq . 

The accretion luminosities (both from disk and hot- 
spot) are turned on in Model 5. However, since the disk 
is truncated at the dust destruction radius which in this 
case is ~ 98.4 r*, most of the disk accretion luminosity is 
lost and the rest of it (Ldisk = 25.29 Lq) is much lower 
than the hotspot luminosity (I/hotspot = 2.45 x 10^ Lq). 
Note that in Robitaille model, part of this missing disk 
luminosity is included in the stellar + hotspot luminosity. 

In Model 6 we adjust the disk to be a geometrically 
thick one {Hq/t^ = 0.1) to see the effects of the height 
of the disk. 

In our fiducial model, the disk is extended to the stel- 
lar radius so that the accretion luminosity can be con- 
served. Inside the dust destruction front (T > 1600 K), 
gas opacities are used. Here we assume that the a-disk 
model is still valid. The scale height, i^o, and radial and 
vertical scaling parameters a and P are calculated self- 
consistently from the assumed opacities and other stellar 
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TABLE 1 

Properties of features included in the model series. 



Models 


Star 


Disk 


Envelope 


Outflow 


1 




no 


52M0, p oc r-^-^ 


no 


2 


8Mq 


no 


52M0, expansion wave 


no 


3 


8Mq 


l/3m*, thin (Ho/r^ = 0.01), 
passive, rd,in = r^uh 


49.333M0, expansion 
wave, rotation 


no 


4 


8Mq 


as above 


only ^ 29M0 left, 
expansion wave, rotation 


yes 


5 


8Mq 


l/3m*, thin, 
rd,in = ^sub, active 


as above 


yes 


6 


8Mq 


l/3m*, active, rd,in = '^sub 
thick (i/o/r* = 0.1) 


as above 


yes 


7 


8Mq 


same as Model 8, 
except rd,in = rsub 


as above 


yes 


8 

(fiducial) 


8Mq 


l/3m*, active, 
ifo/r* = 0.06, rd,in = r* 


as above 


yes 



TABLE 2 

Parameters of the models comparing different column density. 





Model 8 (fiducial) 


Model 81 


Model 8h 


mean surface density Ed (g/cm^) 


1 


0.316 


3.16 


core radius i^core (pc) 


0.057 


0.10 


0.032 


position of expansion wave front rew (pc) 


0.049 


0.088 


0.028 


formation time tf (yr) 


1.29 X 10^ 


3.06 X 10^ 


5.44 X 10^ 


position of sonic point rsp (AU) 


2.57 X 10^ 


4.57 X 10^ 


1.45 X 10^ 


outer radius of disk r^z (AU) 


449.4 


801.4 


253.4 


disk accretion rate rh (Mq /yr) 


2.398 X 10-4 


1.035 X 10-4 


5.667 X 10-4 


stellar radius r* (Rq) 


12.0 


11.3 


5.93 


stellar surface temperature T* (K) 


1.22 X 10^ 


1.25 X 10^ 


1.74 X 104 


stellar + hotspot temperature T^,hotspot (K) 


1.42 X 10"^ 


1.37 X 10^ 


2.62 X 104 


stellar luminosity L* (Lq) 


2.81 X 10^ 


2.82 X 10^ 


2.84 X 10^ 


total accretion luminosity Lace ( ^0) 


4.90 X 10^ 


2.25 X 10^ 


2.36 X 104 


disk scale height at stellar surface if 


0.06 


0.053 


0.078 



and envelope parameters: Hq/v^ = 0.06, a = 1.75, and 
/3 = 1.08. Fig. [2] is a schematic of the structure in the 
fiducial model. The positions of core boundary, expan- 
sion wave front, sonic point, outer radius of the disk, dust 
destruction radius, and the scale height of the disk are 
all marked. Unlike the previous models, in which disk 
material can exist at any height (though the density is 
very low at large z), we truncate the disk at a fixed den- 
sity of 2 X 10~^^ g cm~^ so that the disk atmosphere 
joins smoothly with the infalling envelope. (The input 
density profile can be seen in the lower panel of Fig. [3l) 
We expect that the existence of optically thick gas in the 
innermost region of the disk can make a significant dif- 
ference. In order to investigate this effect, we construct 
Model 7 to be the same as the fiducial model except 
without gas disk {vd^in = ^sub)- And the fiducial model 
is labeled with Model 8. Table [1] lists the differences of 
all these models. 

The column density of the clump where the core is 
embedded in can affect the surface pressure of the core, 
therefore its size, outer radius of the disk, accretion rate 
and also the scale height of the disk. In order to inves- 
tigate this effect, we consider other two values for the 
mean surface density of the clump, Sd = 0.316 and 3.16 
g/cm^. These correspond to a change in column density 
of the core by a factor of 10 and the surface pressure by a 
factor of 100. These two models are labeled with Model 



81 (low surface density) and 8h (high surface density). 
Table [2] lists the differences of these three models. 

3. SIMULATIONS 

We use the Monte Carlo radiation transfer code by 
IWhitnev et al.l (|2003bf ) to perform our calculations. This 
code includes thermal emission, non-isotropic scattering 
and polarization due to scattering from the dust in a 
spherical-polar grid, using the method of Monte Carlo ra- 
diative equilibrium by Bjorkman & Wood (2001). This 
method requires that the opacities are set up at the be- 
ginning of each run and kept unchanged. During the run, 
the temperature in each cell is always being corrected by 
calculating the balance of absorption and emission of new 
photon packets. This becomes a problem when gas opac- 
ities are included because they are highly dependent on 
the temperature and density. So we have to iterate to 
obtain the correct temperature profile for the disk (es- 
pecially for the inner region) and then use it to set up 
opacities for each cell. 

We choose the analytical temperature profile from a- 
disk model in our case as an initial condition. After sev- 
eral iterations, the temperature of the outer disk > 4 AU 
becomes quite stable. However, because of the disconti- 
nuity of the opacities at the transition region between 
gas and dust, the temperature profile in the inner re- 
gion oscillates between two phases. Even though the dif- 
ference can be large between results from two adjacent 
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Fig. 3. — Input temperature and density for our fiducial model. riHe = O.lrin is assumed here. The black contours in the upper panels 
correspond to the dust destruction front (1600K), cells inside which are assigned with gas opacities. The white contours show the photon- 
diffusive region calculated in the fiducial model, a photon emitted inside which will move to the surface through a path with constant r 
and emitted again from the surface. The temperatures of the cells on that path are calculated with grey atmosphere approximation, which 
is why the the dust destruction fronts have these vertical structures. Also, since the frequency of the photon, when it finally leave the 
photon-diffusive region, is calculated from the surface temperature, these vertical structures inside it should not affect the results. 



steps (i.e. Tn{r^O) and Tn-i(r, 6>)), the averaged tem- 
perature profile of two adjacent step should keep similar 
(i.e. log(T^(r,^)) = (log(T„(r,^))+log(T„_i(r,^)))/2). 
So we stop the iteration when ~ ^n-i^ more specif- 
ically, when standard deviation of the distribution of 
(logT^— logT^_;L)/logT^ for the cells inside 4 AU becomes 
consistently < 0.1, which corresponds to a variation of 
^ 25%. Then we use the temperature profile which is 
higher in the midplane to be the input for the final run 
with large number of photons to generate SEDs and im- 
ages. We performed simulations with both input tem- 
perature profiles and saw no significant differences on the 
SEDs and images. We also doubled the photon number 
for the iteration but found no significant dependence of 
the standard deviation on the photon number. It should 
be noted here that 0.1 standard deviation inside 4 AU is 
only an arbitrary standard. The upper panel of Fig. [3] 
shows the input temperature profile. The black contours 
correspond to the dust destruction front (T = 1600 K), 
inside which gas opacities are assigned depending on the 
temperature and density. We can also see that the dust 
destruction front extends to ~ 10 AU in the midplane 
and ^ 3 AU on the surface of the disk, which agrees 
with the estimate of the dust destruction radius of 5.5 
AU in Model 3-7. We note that temperature iteration 
is only used in Model 8, while in other models, the dust 
opacities are assigned only depending on the region and 
the density. Also, this input temperature is only used to 
assign opacities. It is not the initial temperature condi- 
tion for each run. 

The inner region of the disk around the midplane is 
very optically thick and so detailed radiative transfer 
simulation becomes very time-consuming. Thus in the 



code the grey atmosphere approximation is used to de- 
scribe this region (photon-diffusive region). A photon 
generated inside this photon-diffusive region will move to 
the surface of the region by following a path with same 
r, and is then emitted with a frequency calculated based 
on the temperature of the surface cell. The temperatures 
of the cells on that path are calculated accordingly with 
the grey atmosphere approximation. In models except 
the fiducial one (Model 3 - 7), as in the original code, 
a cell is set to be in the photon-diffusive region if the 
optical depth from z = oo to it is larger than 10. In the 
fiducial model, we change this to a more restrictive local 
definition, that the photon-diffusive region is where the 
mean free path is smaller than 0.1 H{r). Note that as 
long as this photon-diffusive region is set small enough, 
it should not affect the final results. 

A 3000 X 1499 grid is used to resolve the r x 9 space. 
In r space, ~ 600 cells are used to resolve the disk with 
a finer grid in inner region, ~ 150 cells are used to cover 
the region inside ~ 6 AU in the fiducial model. In 6 
space, the grid is finer in the disk (especially around the 
midplane) and around the opening angle of the outflow 
cavity O^^qsc- 700 cells are used between 20° above 
and below the disk midplane and ~ 300 cells are used 
between 50° and 53° (<9^,esc = 51°). 

For each run, SEDs at ten inclinations (evenly dis- 
tributed in cosine space) can be produced simultaneously, 
while if the "peeling-off" mode is used, images and SEDs 
with higher signal-to-noise ratio are produced for the par- 
ticular "peeling-off" angle. The "peeling-off" mode is 
very time-consuming. For most of our models, 10^ pho- 
tons are used in one run, and it typically takes several 
days to one week running on a single processor. This 
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Fig. 4. — SEDs for the model series, assuming a distance of 1 kpc. 
A typical inclination of 60° is chosen. The upper panel shows the 
SEDs and the lower panel shows the differences of each model from 
the fiducial one. The black solid curve is for the fiducial model, 
the red thin curves are for Model 1-3 (models without outflow 
cavities), and the blue thick curves are for Models 4-7 (models 
with outflow cavities) The black dotted line is the input stellar 
spectrum (black-body) which is used for all the models. 

number of photon pacl^ets is still not perfect for an im- 
age of the whole core, especially for those at wavelengths 
with low fluxes. Because this code does not enable paral- 
lel computing now, in order to save time, for each model 
we run 10 times with different random seeds simultane- 
ously on different processors, and superposed their re- 
sults, making our final images contain 10^ photons. 

Images are produce at several wavelengths and con- 
volved with filter functions for comparison to observa- 
tions. The code has already included filters such as 
Spitzer IRAC filters at 3.6, 4.5, 5.8 and 8.0 /im, Spitzer 
MIPS filters at 24, 70 and 160 /im, and 2MASS J, H, K 
bands. We also add in the filters of GTC-CanariCam, 
Herschel PACS and SPIRE, and SOFIA FORCAST. We 
will show images both before and after convolution with 
resolution of these particular instruments. 

3.1. Dust Opacity 

We use three dust grain models for different regions: 
(1) the envelope; (2) lower density regions of the disk, and 
(3) the densest re gions in the disk (nn, > 10^^ cm~^, the 
criterion used by iRobitaille et al"] [2006). For our present 
models there is no dust (or gas) in the outflow cavities. 
The default d ust models in the c ode are used without 
any change (Whi tnev et"aIl l2003aV 

The dust model used in the envelope is based on that 
derived by iKim et ahl (|l994l ) for the diffuse interstel- 
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Fig. 5. — The same as Fig. [4l except for an inclination of 30°. 



lar medium (ISM). Th e size distribution is modified by 
I Whitnev et al.l (|2003a[ ) to fit an extinction curve typi- 
cal of the more dense regions of the Taurus molecular 
cloud with Ry = 4. These grains also include a water 
ice mantle covering 5% of the radius. For the lo wer den- 
sity re gions of the disk, the grain model of Coter a et al.l 
(|200ll ) is used. It has grains larger than ISM grains, but 
not as large as the disk midplane grains. 

For the densest regions of the disk, we use the dust 
model with large grains (~1 mm; model 1 in Wood et al.l 
[2QQ2h . which fit the HH 30 disk SED well. Com- 
pared to ISM grains, the larger dust grain model has a 
shallower wavelength-dependent opacity: lower at short 
wavelengths and larger at long wavelengths. 

3.2. Gas Opacity 

For most regions of the disk, dust grains dominate the 
opacity, even if the mass ratio between dust and gas is 
low (0.01 is used in our models). However, in the in- 
nermost region of the disk, where the dust cannot exist 
(T > 1600 K), opacity is dominated by gas. Especially 
when the temperature is high (~ 10^ K), the mean opac- 
ities of the gas can be comparable or higher than dust 
opacities. As discussed in Section 12. 2[ most of the disk 
accretion luminosity is from this innermost region. Thus, 
gas opacity should be included in a realistic model. 

In order to simulate the frequency of each photon 
packet emitted from the gas-dominated region correctly, 
not only Rosseland or Planck mean opacities, but also 
the frequency-dependent opacities are needed. Besides, 
the gas opacity is highly dependent on the temperature 
and the density, which make our problem much more 
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Fig. 6. — Output temperature distribution of Model 7 and 8, showing the difference the gaseous innermost disk makes. The black 
contours are dust destruction front (T=1600K). 



complicated. Since our aim is only to simulate the disk 
emission correctly, especially in near and mid-IR, rather 
than to simulate the details of line profiles, we smooth 
the monochromatic opacities for simplicity and smaller 
memory demand. Also, we assign opacities to a grid in 
temperature and density, rather than to interpolate to 
obtain opacities at exact temperature and density. 

For temperatures higher than ^ 3000 K, we adopt 
gas opacities from OP project (ISeaton et ahl 11994 
iBadnell et all l2005l and ISeatonI l2005f ). They provide 
monochromatic gas opacities of hydrogen, helium and 
other 15 elements, in large ranges of temperature (~ 3000 
K to 10^ K, we only use those of temperature up to 
~ 10^ K in our present model since the maximum tem- 
perature we can reach in the models is less than this) 
and density. The opacities are mainly due to the line ab- 
sorption, ionization, H~ bound-free absorption, electron- 
hydrogen/ helium free- free absorption. For scattering, 
Tho mpson scattering with a collective effect is consid- 
ered (|Boercker|[T987n . For regions with T < 3000 K, 
the opacity model of [Alexander fc FergusonI (1994) is 
adopted. However, they do not provide monochromatic 
opacities for a range of temperatures and densities. So 
we use the monochron i atic o pacities shown in Fig. 4 of 
[Alexander fc FergusonI (|l994l ) which is for a temperature 

of 2000 K and a density of 8 x 10~^^ g/cm^, and rescale 
it for other temperatures and densities based on their 
Rosseland mean opacities. Fortunately, the temperature 
range for this model is quite narrow (1600K - 3000K). 
Alexander and Ferguson's model includes both gas and 
molecules. At ~ 2000K, the total opacity is dominated 
by molecules, especially H2O and TiO. Atomic lines and 
CO lines are important at some wavelengths. Other con- 
tinuous sources include H~ absorption and Rayleigh scat- 



tering of H and H2. 

In this way we construct an opacity grid of temperature 
and density. The temperature varies from 1600 K to 10^ 
K with an interval of 0.1 in logarithmic space, and the 
density varies from ~ 10~^ to ~ 10~^^ g/cm^ with an 
interval of 0.5 in logarithmic space. For each T and p, a 
monochromatic opacity file is assigned, so that we have 
totally 148 gas opacity files in our present model. At the 
beginning of each run, a temperature profile is read in, 
and in each cell the gas opacity file with closest T and 
p to the read-in values is chosen. The opacities are not 
changed during the calculation. 

Since here we are not concerned with line absorption 
and emission, it is better to smooth the monochromatic 
opacities to save computing memory. In the code, three 
important values are calculated: 1) Rosseland mean 
opacity 

_1_ ^ j^{dB,ldT)/nUv)Av 

«R /o°°((9i?,/(9r)di. ' ^ ' 

used to determine the photon-diffusive region inside the 
disk; 2) Planck mean opacity 

Jo ^^^^ 

used to calculate the energy equilibrium in a cell; and 3) 



(32) 



/o /^p dT 
where Py is the probability that a photon packet is 
emitted from a cell with a frequency between and v 
(|Biorkmanl Il997[ ) . For tip and P^, it is best to smooth 
opacities with linear averaging, thus giving better esti- 
mates of equilibrium temperatures and photon frequen- 
cies. However, this method tends to increase kr by 
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Fig. 7. — SEDs of the fiducial model (Model 8) at 4 inclinations, 
assuming a distance of 1 kpc (F^ in upper panel and uFu in lower 
panel). The thicker lines are total fiuxes and the thiner red lines 
are only the scattered light. The lower dotted blue line is the input 
stellar spectrum (black-body). The upper dotted line is the total 
luminosity of star and hot-spot emitted as a black-body. 



~10% or more. Averaging \ogK reduces the importance 
of line absorption and yields more accurate Rosseland 
mean opacities. For our problem, the precise location 
of the boundary of the photon-diffusive region is is not 
very important, so we smooth opacities with linear aver- 
ages. The original gas opacity files contain 10^ frequen- 
cies between hv/kT = 0.001 and 20. We smooth them 
by averaging 50 adjacent frequencies. 

4. RESULTS 

4.1. SEDs 

4.1.1. SEDs of the Model Series 

Figures 3] and [5] show the SEDs of the model series at 
inclination angle between our line of sight and the pro- 
tostar rotation axis of 60° (i.e. a more equatorial view) 
and 30° (i.e. a more polar view) respectively. In each 
figure, the upper panel shows the total fluxes while in 
the lower panel these SEDs are compared to the fiducial 
model (Model 8). The = 60° line of sight goes through 
the envelope material, in which case the photons from the 
protostar and the disk are all reprocessed by the dust be- 
fore they can escape the core. The = 30° line of sight 
from the central star passes through the outflow cavity 
(for Model 4 - 8), in which case we can see the stellar 
black-body spectra in the short wavelength region of the 
SEDs. Some important features in the SEDs include the 



water ice feature at 3 /im and silicate features at 10 /im 
and 18 /im. The ice feature is only present for the higher 
inclination view for which the lines of sight pass through 
the envelope material, which uses a dust model with ice 
mantles. 

Models 1 to 3 show very similar SEDs, where we do 
not see any radiation at short wavelengths. The 10 /im 
silicate features are all very deep. The occurrence of the 
expansion wave (Model 2) and rotation/disk (Model 3) 
shift the far-IR peak a little to shorter wavelengths, and 
increases the mid-IR emission, making the 18 /im sili- 
cate feature deeper. This is because the expansion wave 
and the Ulrich solution decrease the density of the inner 
region of the core, and thus reduce the extinction. For 
Model 3, which is not spherically symmetric, the SEDs 
do not show much difference between the different in- 
clinations. It should be noted that in Model 3 rotation 
is only considered inside the sonic point. In a more re- 
alistic solution, the material in the outer region of the 
envelope would also be redistribute d by the effect of the 
rotation (like the solution of Terebe v et al.l ([l984) for an 
isothermal core) so that one might see larger differences. 

The outflow cavities change the shape of SEDs signif- 
icantly. With outflow cavities (Model 4 to 8), the SEDs 
show a large excess at wavelengths shorter than 10 /im. 
The position and height of the far-IR peak and the 20 
/im ~ 70 /im slopes are affected by the cavity as well. 
Especially for a low inclination, the star and disk can be 
seen directly. The fluxes at wavelengths shorter than 10 
/im are larger than those at a high inclination by about 
two orders of magnitude. Note that the short wavelength 
emission seen in the = 60° view is essentially all due to 
scattered emission from the outflow cavity walls. 

Compared to a passive disk, an active disk with ac- 
cretion luminosity increases the fluxes at all wavelengths 
without many changes in the shape of the SEDs (compar- 
ing Model 4 and 5). The accretion luminosity in Model 
5 is mainly due to the hot-spot, while most of the disk 
accretion luminosity is lost here because of the absence of 
the innermost disk. The total energy is conserved only in 
Model 8, which the disk is extended to the stellar surface. 

The effect of the thickness of the disk is distinct on the 
SEDs (comparing Model 5 and 6). A thicker disk tends 
to obscure more photons at high inclinations and emit 
or scatter them to low inclination directions - i.e. more 
flux escapes via the outflow cavities. Therefore, with a 
thicker disk. Model 6 at ^ = 30° shows a rise between 
1 and 10 /im which even smooths out the far-IR peak, 
while at ^ = 60° the SED shows a decrease in near-IR 
and shorter wavelengths, and a deeper silicate features. 

Model 8 shows how the SED changes by including the 
flux from the innermost disk region. Compared to pre- 
vious models, the optical and near-IR emission is sig- 
nificantly increased, in both high and low inclinations. 
Model 7 has exactly the same geometry and density setup 
as Model 8, except that it has a disk truncated at rgub- 
The opacities are chosen depending on the input tem- 
perature and density in Model 8 while in Model 7 they 
only depend on the density, therefore, even outside rgub 
the opacity setup of these two models can be different. 
Because of the truncated disk, the majority of the disk 
accretion luminosity is lost in Model 7. This missing 
luminosity shows up in Model 8 mainly as optical radi- 
ation. However, the near-IR flux in Model 8 is not so 
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Fig. 8. SEDs of Model 8 (Eci= 1 g/cm^), Model 81 (Eci= 0.316 g/cm^) and Model 81i (Eci= 3.16 g/cm^) at inclinations of 60° (left 
panels) and 30° (right panels), assuming a distance of 1 kpc. For each inclination, the upper panel shows the total fluxes (black), scattered 
fluxes (red) and stellar input (blue), while the lower panel shows fluxes from the disk only. 



bright as Model 7. The far-IR SEDs of these two models 
do not show much difference. 

Fig. [6]shows the final output temperature distributions 
of these two models. The temperature of the disk reaches 
~ 10^ K at the innermost region in Model 8, while the 
highest temperature in Model 7 is only ~ 3000 K. This 
explains the optical excess in the SED of Model 8. How- 
ever, the existence of the innermost disk shields the fiux 
from the star, therefore in Model 7 the temperature at 
the surface of the disk and the base of the outfiow cavity 
becomes higher, leading to the higher near-IR fiux on the 
SED. For the models presented here, we do not have any 
material in the outfiow cavity. If any dust grain exists 
there, the optical emission in Model 8 will be suppressed 
and the disk luminosity may appear as more near- and 
mid-IR radiation. This will be examined in a future pa- 
per. 

4.1.2. SEDs of the Fiducial Model 

Fig. [3 shows the SEDs of our fiducial model (Model 
8) at four inclinations. Even after we smoothed the 
monochromatic opacity curves of gas, the original out- 
put SEDs still show some significant emission features, 
mainly Ha. Since the line profile is not exactly correct 
in our models and it is not the interest of our present 
study, we subtract the Ha feature and smooth out other 
line features with wavelengths < 2 /im. The energy con- 
tained in the Ha line is typically < few%. 

The SEDs at 4 inclinations are shown here, from a view 
along the rotation/outfiow axis to that through the disk 
plane. The inclination of the viewing angle changes the 



observed fiux at short wavelengths significantly. It also 
affects the height and position of the peak in the far-IR 
region, and the mid-IR spectral slope. The SED at wave- 
lengths longer than ~100 /im is not affected by the incli- 
nation. Inclinations of 30° and 0° have very similar SEDs 
(0° inclination contains 8% more energy than 30° incli- 
nation). Recall, the outfiow cavity has an opening angle 
of 51°, which means we can directly see the star in both 
these cases. The stellar spectrum and the black-body 
spectrum containing both stellar and hot-spot luminos- 
ity are shown by the dotted lines. Their difference shows 
the luminosity from the hot-spot region. From high in- 
clinations to low inclinations, because of the change of 
the optical depth, the silicate feature change from a big 
absorption feature to a weak emission feature. 

The dashed lines show the SEDs of only scattered light. 
At high inclination angles, the observed light at short- 
wavelengths has always been scattered, while at low in- 
clinations, we can see stellar radiation and thermal emis- 
sion from the disk. In the far-IR, the fiux is dominated 
by the thermal emission of the envelope. Such significant 
difference between low inclinations and high inclinations 
is partly because we have an empty outfiow cavity. 

4.1.3. Effect of Different Mass Surface Density 

Model 8, 81 and 8h compare the effect of different mean 
surface densities of clump in which the core is embed- 
ded, which affects the surface pressure of the core, and 
therefore the size of the core, the accretion rate, the disk 
structure and the protostar evolution. The size of the 
core, the radius of the expansion wave front, the position 
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Fig. 9. — Images for the fiducial model at inclination of 60°, assuming a distance of 1 kpc, at different wavelengths. Each image has 
149 X 149 pixels, and field of view of 30" x 30". Each image is normalized to its maximum surface brightness which is labeled in the 
bottom-left corner. The values on the top-right corners are the total fiuxes. The blue contour is 1% of the maximum surface brightness 



of the sonic point and the disk radius are ah propor- 

— 1/2 

tional to , so with higher clump surface density, 

the core is more compact and accordingly the accretion 
rate is higher. The scale hight of the disk is also larger in 
the high Ed case. The stellar luminosities of these three 
models are similar but Model 8h has a much bluer stel- 
lar spectrum, because of the higher hot-spot accretion 
luminosity. Some important parameters of these three 
models are listed in Table. [2] 

Fig. [8] compares SEDs of these three models at both 
60° and 30°. As discussed above, higher surface den- 
sity leads to higher bolometric luminosities, which can be 
seen from the SEDs: At inclination of 30°, with a higher 
surface density, the flux is higher at all wavelengths; And 
at 60° inclination we can see the same effect in opti- 
cal, near-IR and far-IR emissions. However, in mid-IR, 
Model 81 shows a rise of the flux and a very significant 



10 /im silicate feature, while the other two models do 
not. To explain this, we also show the disk flux in the 
lower panels. Here, disk flux contains photons which have 
their last emission in the disk, and then either escape the 
core directly or are scattered before they reach the ob- 
server. At the inclination of 60°, the disk is blocked by 
the envelope. The short- wavelength fluxes from the disk 
should all have been scattered. They keep the trend that 
the model with higher surface density has higher fluxes. 
However, in mid- and far-IR, the fluxes should have suf- 
fered the extinction of the envelope, making the model 
with higher surface density have lower fluxes because of 
the higher extinction. Thus, even though Model 8 and 
8h have very strong silicate features in their disk SEDs, 
they are buried in the envelope fluxes in the total SEDs, 
while Model 81 shows the high disk flux level in mid-IR 
and a significant silicate absorption feature because of its 
lower extinction and lower envelope flux. 
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Fig. 10. — Same as Fig. [O] except convolved with the resolution of each particular instrument, labeled in the top-right corner and shown 
as red circles of radius equal to the HWHM (or line segments with lengths equal to the HWHMs for Herschel SPIRE) on the bottom right 
corners. The maximum surface brightness (which is dependent on the resolution) is labeled in the bottom-left corner. The contours for 
most of the images have intervals of half an order of magnitude from the maximum brightness. Because some images (1.23, 10.3 and 20 
)Lim) are relatively noisy at faint fluxes, only flve contours with intervals of half an order of magnitude are shown here. For images at 250 
and 500 /im, dashed contours have intervals of 0.1 dex. 



4.2. Images 
4.2.1. Images of the Fiducial Model 

Fig. [9l [TOl [11] and [12] show the images for the fiducial 
model (Model 8) at inclinations of 60° and 30°, assuming 
a distance of 1 kpc and no foreground extinction. The in- 
struments filters chosen here are: 2MASS J, H, K bands, 
Spitzer IRAC 3.5 /im, 4.5 /im, 5.8 /im and 8 /im, MIPS 24 
/im, GTC-CanariCam N band and Q wide band, SOFIA 
FORCAST 37.1 /im, Herschel PACS 70 /im, 100 /im, 160 
/im, and SPIRE 250 /im and 500 /im. We show both re- 
solved images and those after convolving with PSFs of 
the particular instruments. Images showed here are all 
normalized to their maximum surface brightnesses which 
are labeled on each image. 10^ photon packets are used. 



However they are still not enough to produce a smooth 
image for the narrow filters or at the wavelengths with 
low fluxes. Besides, the simulation grid also contributes 
to this problem, like the strip patterns in the images at 
some wavelengths (e.g. 20 /im). Around the opening an- 
gle of the outflow cavity, the grid with very small change 
of the polar angle intersects with the cavity wall at quite 
different radii. Even we have made the grid much finer 
in polar angle at the region around the cavity wall, these 
patterns can still be seen, especially at outer regions of 
the envelope. A finer grid would demand more memory 
and computing time. 

On the images at inclination of 60° before convolu- 
tion with the instrument PSF, the most significant fea- 
tures are the outflow cavities. They can be seen in any 
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Fig. 11. — Same as Fig. O except at inclination of 30° 



wavelength, though they become not so obvious in far-IR 
wavelengths where the thermal emission of the envelope 
dominates. At this inclination, the line of sight intersects 
with the envelope, thus in near-IR the only emission we 
can see is scattered by the cavity wall and escapes from 
the opening region, especially from the side facing us. 
Deeper regions appear as the wavelength increases. The 
central protostar and disk begins to show up in mid-IR 
images. The dark lane around it shows the size of the 
disk. In the mid-IR, the cavity walls do minates the emis- 
sions, as has been discussed bv lDe Bu izer (2006). Both 
sides can be clearly seen and the brightest region is the 
base of the cavity. The emission of the envelope takes 
over at longer wavelengths, making the image symmetric 
and the outflow cavity begins to fade. The central star 
and disk can still be seen as the brightest region in those 
wavelengths. 

At the inclination of 30°, the central object can always 



be seen through the empty outflow cavity. At shorter 
wavelengths, only the side facing us of the cavity is signif- 
icant and the opposite side is very dim. Far-IR emission 
is dominated by the envelope and it is hard to tell the 
features of the cavities. It should be noted here, espe- 
cially when comparing the 30° and 60° images, that the 
images are all scaled to their maximum surface bright- 
ness, which is generally much greater for the 30° viewing 
angles. 

After convolving with PSF of particular instruments, 
in the far-IR, the contours become very symmetric and 
we cannot tell the inclination or the opposite outflow 
cavities. It is possible to see the opposite outflow in MIPS 
24 /im and FORCAST 37.1 /im at high inclinations, if the 
S/N is large enough. The two sides of the outflow cavities 
are clear in IRAC images and in near-IR images. GTC- 
CanariCam has very high resolution so it may enable us 
to see the central region in mid-IR. 
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Fig. 12. — Same as Fig. 1101 except at inclination of 30° and that for images at 1.23, 10.3 and 20 /im, only five contours with intervals of 
an order of magnitude are shown. 



Fig. [13] shows the flux profile along the system axis 
at inclinations of both 60° and 30°. Each profile is nor- 
malized by the mean flux of a very narrow strip along 
the axis across the whole core. Thus, the general shape 
and level of the profiles are independent on the resolu- 
tion of the image. But the curve would be smoothed out 
if a larger PSF-FWHM is used. At inclination of 60°, 
as we discussed above, at short wavelengths, the central 
region is totally blocked by the envelope, most of the 
emission comes from the outflow cavity. In the mid-IR, 
the maximum surface brightness comes from the base 
of the outflow cavity rather than central star and disk, 
on the side the outflow cavity is opposite to us, the flux 
drops very fast by almost four orders of magnitude, while 
on the side the outflow is facing us the flux drops much 
more gradually. In the far-IR, the profile is symmetric on 
both sides, and central region has the maximum surface 
brightness. At lower inclination, the maximum flux al- 



ways comes from the center. In the far-IR, the profile is 
very similar to the case in higher inclination. At shorter 
wavelengths, the contrast between outer regions and the 
center point is much higher than the previous case, so 
that the profile drops very fast. At 1.66 fim and 10.3 fim 
we can only see the flux due to the outflow cavity facing 
us, while the other side only shows up at 20 jam. Because 
in our models it is empty outside the core, the profiles 
all cut off at the core radius, which is not true in reality 
since clump material can extend to larger radii. 

4.2.2. Effect of Different Surface Density 

Fig. [HI compares the images of Model 8, 81 and 8h at 
6 wavelengths. The size of the images are all 50"x50". 
We can see that the core is smaller when the surface 
pressure is higher. The surface brightness is dependent 
on the resolution, so here we convolved all images with 
same PSF with FWHM of 0.5". The images are scaled 
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Fig. 13. — Flux distribution along the projected outflow axis. 
Fluxes are normalized by the mean value of a narrow strip across 
the core along the axis. Different curves show different wave- 
lengths. The upper panel is at inclination of 60° and the lower 
panel is at 30°. 

to the maximum surface brightness of Model 8h at each 
wavelengths, so the brightnesses of the three images at 
each wavelength can be directly compared. 

The total fluxes and the maximum surface brightness 
in images of Model 8h are generally higher than those 
of the other two models. This is natural since the to- 
tal luminosity is much higher in this model. The optical 
depth is larger when is higher. In the near-IR, we can 
see the base of the outflow cavity facing us and most of 
the opposite side in Model 81, while we can only see the 
light from the opening region of the cavities in Model 
8h. The central region shows up in Model 8h only in 
far-IR, while it can be seen at mid-IR in Model 8 and 
81. With lower surface density, the contrast between two 
sides of the outflow cavities becomes smaller. In Model 
8h the side towards us is much brighter than the other 
side, even at 70 /im, while the other two models show al- 
most symmetric images. We can also see that the mid-IR 
emission is dominated by the cavity walls, especially in 
the cases with higher extinction, as opposed to that with 
lower extinction the emissions are more concentrated to 
the central region. 

5. DISCUSSIONS AND CONCLUSIONS 

We have constructed a model for individual massive 
star formation, with a 60 Mq initial core that forms an 
8M0 protostar. We included the inside-out expansion 
wave in the core, free-fall rotating collapse in the inner 
region, an accretion disk of 1/3 the mass of the star, 
and bipolar outflow cavities with large opening angle, 
parameters of which are all calculated self-consistently. 
For the first time, we considered an optically thick inner 



disk with gas opacities which were assigned depending on 
the pre-calculated temperature profile. This inner disk 
enabled us to calculate the emission from the active disk 
with correct total luminosity and spectrum. 

Compared to the Robitaille model, our parameter 
space covers higher accretion rates, higher disk mass, 
denser envelopes and larger outflow cavities. Also, in the 
Robitaille model the disk accretion luminosity is much 
lower than that in our model since the disk is truncated at 
the dust destruction radius, while the hot-spot luminos- 
ity on the stellar surface is set to be Gm^rh{l/r^ — 1/rco) 
where Tco is the magnetic co-rotation radius and set to 
be 5r*, making the hot-spot luminosity is 4/5 of the to- 
tal accretion luminosity Lace- So energy of ^ 0.2 Lace 
is lost. In our model, the total accretion luminosity is 
divided equally into disk accretion luminosity (emitted 
from the disk extended to stellar radius with gas and dust 
opacities) and hot-spot luminosity (added on the stellar 
spectrum). Also, the Ulrich solution is used for the whole 
envelope in Robitaille model, which means the whole en- 
velope is assumed to be free-falling and the infalling rate 
is constant at all radii. As discussed in Section [2. 1.31 the 
core undergoes free-fall only in the inner region. There- 
fore, for a given core mass and accretion rate, our model 
has a more compact core with higher extinction, causing 
the far-IR peak of the SED to be at lower fluxes and at 
longer wavelengths. 

We have presented SEDs of the model series, the fidu- 
cial model, and the models with higher and lower surface 
pressure, at typical inclinations. We also have presented 
images for the fiducial model at JHK bands, IRAC and 
MIPS bands, GTC-CanariCam bands, SOFIA bands and 
Herschel bands, both resolved and convolved with the 
resolutions of each instrument. The main conclusions 
can be summarized as: 

1. Outflow cavities affect the SEDs significantly and 
cause a large difference between low and high inclinations 
of our viewing angle. However, the present modeling as- 
sumes these cavities are optically thin, which may not 
be valid, especially at the shorter wavelengths. This is- 
sue will be investigated in a follow-up paper. The height 
of the disk also affects the SEDs. With a thicker disk, 
the near- and mid-IR fluxes at low inclinations become 
higher, while at high inclinations it suppresses the fluxes 
(especially at short wavelengths) and creates deep sili- 
cate features. Also, the density distribution in the core 
(especially the inner region) can affect the mid-IR flux 
levels, the silicate features, and the far-IR peaks. 

2. The temperature of the innermost region of the disk 
can reach ~ 10^ K. The disk becomes optically thick in 
such conditions even if no dust can exist there. SEDs 
show the rise of optical emission due to this hot inner 
disk. This optically thick inner disk also shields flux 
from the protostar, leading to lower temperatures on the 
surface of the disk and the base of the cavity wall, and 
therefore lower fluxes in near-IR and mid-IR part of the 
SEDs. 

3. The SEDs of the fiducial model show that the in- 
clination can affect SEDs at wavelengths shorter than 
100 /im, including the far-IR peak. The mid-IR spec- 
tral slope changes significantly with inclination. The sil- 
icate feature changes from a deep absorption feature to 
an emission feature from high inclinations to low inclina- 
tions, due to the change of the optical depth. 




Fig. 14. — Images for Model 8, Model 81 and Model 8h at the inclination of 60°, assuming a distance of 1 kpc. The size of the image is 
50" X 50". The resolution here is set to be 0.5", so the surface brightnesses can be compared. Each images are scaled to the maximum 
brightness of Model 8h at that band, but their own maximum brightness is also labelled in the bottom-left corners. The total flux in each 
image is also labelled in the top-right corner. The blue contour is 1% of the maximum surface brightness in that image. 



4. With higher surface pressure, the core becomes more 
compact and the accretion rate and luminosity become 
higher, leading to higher fluxes at all wavelengths (ex- 
cept in the mid-IR for high inclination cases which suffer 
higher extinction). High extinction can also cause the 
mid-IR flux to be dominated by the envelope, and thus 
hide the silicate absorption feature. Thus, only in the 
model with low extinction does the silicate absorption 
feature show up. 

5. Outflow cavities are the most significant features 
on the images, except at wavelengths longer than 70 /im. 
At inclinations of 60°, from short wavelengths to long 
wavelengths, the brightest point moves from the outer 
region of the cavity to the base of the cavity wall, and 
to the center of the core in far-IR, while at inclination 
of 30°, it is always the central region. At inclination 
of 60°, the opposite outflow cavity can be seen n the 
mid-IR if fluxes ~ 10^ — 10^ times fainter than the peak 
can be probed. It is very difficult to see the opposite 
cavity at an inclination of 30°. GTC-CanariCam (and 
other 10 m class telescopes with mid-IR cameras) has 
very high angular resolution so it may enable us to resolve 
the central disk system. The flux distribution along the 
outflow axis can help constrain model assumptions and 
inclination angles and so will useful to measure. In the 
mid-IR, the cavity walls seem to dominate the emission. 



but for lower density cores with lower extinction, the 
central region becomes brighter. The contrast between 
the two sides of the outflows in mid- and far-IR increases 
as the extinction becomes higher. 

The model we have presented will be improved in fu- 
ture papers by a detailed consideration of the effect of 
including gas and dust in the outflow cavity. Currently 
models for the density distribution here are quite uncer- 
tain, which is why we defer this study to a future date. 
An evolutionary sequence of protostellar models with and 
without outflow opacity will then be presented. To gauge 
the degree of inhomogeneity we will consider the re- 
sults of numerical sim ulations of core accretion and out- 
flow interaction (e.g. iKrumholz et al.l [20071 IStaff et all 
2010), and then incorporate these inhomogeneities into 
the models in a parametrized form. 
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